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Using high-resolution SPH numerical simulations, we investigate the effects 
of gas on the inspiral and merger of a massive black hole binary. This study is 
+3 • motivated by the very massive nuclear gas disks observed in the central regions 

of merging galaxies. Here we present results that expand on the treatment in 
a previous work (Escala, Larson, Coppi & Mardones 2004; henceforth Paper I), 
by studying more realistic models in which the gas is in a disk with significant 
dumpiness. We run a variety of models, ranging from simulations with a rel- 
atively smooth gas disk to cases in which the gas has a more clumpy spatial 
distribution. We also vary the inclination angle between the plane of the binary 
and the plane of the disk, and the mass ratio between the MBHs and the gaseous 
disk. We find that as in Paper I, in the early evolution of the system the bi- 
nary separation diminishes due to gravitational drag, and in the later stages the 
medium responds by forming an ellipsoidal density enhancement whose axis lags 
behind the binary axis; this offset produces a torque on the binary that causes 
continuing loss of angular momentum and is able to reduce the separation to 
distances where gravitational radiation is efficient. The main difference is that 
between these two regimes we now find a new transition regime that was not 
apparent in Paper I, in which the evolution is temporarily slowed down when 
neither of these mechanisms is fully effective. In the variety of simulations that 
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we perform, we find that the coalescence timescale for the MBH binary varies 
between 5 x 10 6 yr and 2.5 x 10 7 yr. For MBHs that satisfy the observed 'm — cr c ' 
relation, we predict that in a merger of galaxies that have at least 1% of their 
total mass in gas, the MBHs will coalesce soon after the galaxies merge. We also 
predict that if the MBHs depart considerably from the 'm — <7 C ' relation, strong 
tidal and/or resonant forces from the MBH binary can create a circumbinary 
gap in the disk that stalls the coalescence, but this gap formation can act as a 
self-regulatory mechanism on MBH growth that can help to explain the existence 
of the 'm — «j c ' relation. Our work thus supports scenarios of massive black hole 
evolution and growth in which hierarchical merging plays an important role. The 
final coalescence of the black holes leads to gravitational radiation emission that 
would be detectable out to high redshift by LISA. 

Subject headings: Black Hole Physics: binaries, hydrodynamics - Cosmology: 
theory - Galaxies: evolution, nuclei - Quasars: general. 



1. Introduction 

The possibility of massive black hole (MBH) mergers follows from two widely accepted 
facts: that galaxies merge and that every galaxy with a significant bulge hosts an MBH at its 
center (Richstone et al. 1998). This problem was first considered by Begelman, Blandford, 
& Rees (1980) in a study of the long-term evolution of a black hole binary at the center of 
a dense stellar system. Initially, dynamical friction brings the two black holes toward the 
center of the system. The resulting binary MBH continues to shrink via 3-body interactions 
with the surrounding stars, but 3-body interactions tend to eject stars from the central 
region, causing the merger eventually to stall, unless some additional mechanism is able to 
extract angular momentum from the MBH binary. 

The possible additional mechanism that extracts angular momentum from the binary is 
very likely to be gas dynamical in origin. Observations of gas-rich interacting galaxies such 
as the 'Ultraluminous Infrared Galaxies' (ULIRGs) indicate that large amounts of gas are 
present in the central regions of merging galaxies (Sanders & Mirabel 1996). This massive 
nuclear gas concentration has an important influence on the evolution of any central MBH 
binary that forms following a galaxy merger. Since the gas is strongly dissipative, unlike the 
stars, it is expected to remain concentrated near the center and thus to play a continuing 
role in driving the evolution of a central binary MBH. Our aim is to study numerically the 
role of massive nuclear gas clouds in driving the evolution of a binary MBH. 
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In a previous paper (Escala, Larson, Coppi & Mardones 2004; henceforth Paper I) we 
studied numerically the role of a massive spherical gas cloud in driving the evolution of a 
binary MBH. We followed the evolution of the binary through many orbits and close to the 
point where gravitational radiation becomes important. In Paper I, we presented results 
for a relatively simple idealized case in which the gas is assumed to be supported by a 
high virial temperature, so that the gas retains a nearly spherical and relatively smooth 
distribution. We found that in the early evolution of the binary, the separation decreases 
due to the gravitational drag exerted by the background gas. In the later stages, when the 
binary dominates the gravitational potential in its vicinity, the medium responds by forming 
an ellipsoidal density enhancement whose axis lags behind the binary axis, and this offset 
produces a torque on the binary that causes continuing loss of angular momentum and is 
able to reduce the binary separation to distances where gravitational radiation is efficient. 
Assuming typical parameters from observations of Ultra Luminous Infrared Galaxies, we 
predicted that a black hole binary will merge within 10 7 yrs. 

The model presented in Paper I is relatively idealized and, as we describe in §2, obser- 
vations suggest that the gas in merging systems is in a rotating nuclear disk. In this paper, 
we study numerically the role of a massive gas disk, like those seen in the central regions 
of ULIRGs, in driving the evolution of a binary MBH. We run a variety of models, ranging 
from simulations with a relatively smooth gas disk to cases in which the gas has a more 
clumpy spatial distribution. We follow the evolution of the binary through many orbits and 
close to the point where gravitational radiation becomes important. 

We start with a review of the properties of the inner regions of merging galaxies in 
§2. We continue with a description of the assumed initial conditions and the model setup 
in §3. In §4 we study the evolution of a binary MBH in a massive gas disk, exploring the 
parameters that may be relevant for the coalescence of the binary MBH. In §5 we study the 
final coalescence of the binary MBH using higher resolution simulations. In §6 we study 
the criteria for opening a circumbinary gap and the implications for the 'm — o" c ' relation. 
Finally, our conclusions are presented in §7. 



2. Gas in the Nuclei of Merging Galaxies 

Observational and theoretical work both indicate that large amounts of gas can be 
present in the central regions of interacting galaxies, and that this gas can be a domi- 
nant component of these regions. Numerical simulations show that in a merger of galaxies 
containing gas, much of the gas can be driven to the center by gravitational torques that 
remove angular momentum from the shocked gas (Barnes & Hernquist 1992, 1996; Mihos 
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& Hernquist 1996); as a result, more than 60% of the gas originally present in the merging 
galaxies can end up in a very massive nuclear disk with a radius of several hundred parsecs 
(Barnes 2002). Observations of gas-rich interacting galaxies such as the 'Ultr aluminous In- 
frared Galaxies' (ULIRGs) confirm that their central regions often contain massive and dense 
clouds of molecular and atomic gas whose masses are comparable to the total gas content of 
a large gas-rich galaxy (Sanders & Mirabel 1996). 

Downes & Solomon (1998, hereafter DS), using CO interferometer data, show that the 
molecular gas in ULIRGs, at subarcsecond resolution, is in rotating nuclear disks. These 
molecular disks are highly turbulent, with typical velocity dispersions of a — lOOkms -1 . DS 
modeled the CO luminosity as coming mostly from a relatively low-density phase with a gas 
kinetic temperature in the range ~ 65 — 100K. This 'low-density' gas is in a smooth medium 
and not in self-gravitating clouds, and has an average density of 10 3 cm~ 3 . DS estimate that 
only ~ 10% of the gas is in highly opaque and dense (10 5 cm -3 ) self-gravitating clouds that 
are stable against tidal forces, although the exact fraction is still unclear. The derived total 
gas mass is high, typically 5 x 1O 9 M , which is similar to the mass of molecular clouds in a 
large, gas-rich spiral galaxy. However, the total mass in molecular gas is still small compared 
to the dynamical mass derived from the rotation curve, the ratio of gas mass to enclosed 
dynamical mass being about M gas /M dyn = 1/6. In addition to these multi-phase nuclear 
disks, DS also found exceptionally massive clumps with ongoing extreme starburst activity 
having characteristic sizes of only lOOpc, with about 10 9 M Q of gas and an IR luminosity of 
3 x 1O 11 L from recently formed OB stars. 

Since ULIGRs are characterized by an ongoing starburst, feedback effects from star 
formation may play a major role in determining the physics of the ISM. The main feedback 
effects on the gas arise from massive stars and include stellar winds, supernova explosions, 
and the effect of radiation pressure. Recently, Wada and collaborators have made some 
attempts to model this complex environment including all the relevant feedback effects (Wada 
& Norman 2001, Wada 2001). They found a multi-phase ISM in which substructure is 
continuously formed and destroyed, but it reaches a quasi-steady state which is characterized 
by a network of cold (T < 100 K) dense clumps and filaments embedded in a hot (T > 10 6 
K) diffuse medium. 



3. Model Setup 

Having reviewed the observations of molecular gas in ULIGRs, we now describe our 
approach to modeling the complex gaseous environment in the inner regions of merging 
galaxies. 
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In our model, the gas is distributed in a rotating disk with a Mestel surface density 
profile (Mestel 1963) 

xw*» = ^ ID 

and it is initially distributed uniformly between z=±0.5. In our simulations, we use the 
following units : [Mass] = 5 x 1O 9 M , [Time] = 2.5 xl0 5 yr, [Velocity] = 156.46 Km s" 1 and 
[Distance] = 40 pc. In these units, the gravitational constant G is 22.0511, the disk radius 
is 10, and the total gas mass is 1. 

In addition to the gaseous disk, in this paper we include a stellar bulge that stabilizes 
the disk and plays an important role in the evolution of the MBH binary when the binary is 
out of the plane of the disk. The bulge is modeled using a Plummer Law (Plummer 1911) 

where a=200 pc is the core radius and Mb u i ge is the total mass of the bulge. The mass of 
the bulge within r=10 is 5 in the just mentioned units. This choice was made to satisfy the 
observed rotation curves (Downes & Solomon 1998) which imply that the total gaseous mass 
is 1/6 of the dynamical mass. We model the bulge with a collection of 100,000 collisionless 
particles, with a gravitational softening length of 4 pc. For the gaseous disk, we employ 
235,331 SPH particles, and the softening length is also 4pc. 

We first considered a cold massive disk (T=200K, M^isk = 5 x 1O 9 M ) with an isother- 
mal equation of state. This temperature is much smaller than the virial temperature, and 
therefore the cold disk fragments into clumps. These clumps undergo a runaway collapse 
because of the isothermal equation of state, and the collapse stops only at the gravitational 
softening length. The result is the formation of clumps with average densities many orders 
of magnitude higher than the observed ones (Downes & Solomon 1998), and therefore the 
simulated clumps are physically unrealistic. In this unrealistic model, the MBHs interact 
violently and chaotically with the clumps, but in most cases quickly end up at the center 
where they continue to lose angular momentum by interaction with a central massive clump, 
qualitatively as in Paper I. 

In reality, in the inner regions of ULIRGs, feedback processes from starburst activity 
may tend to halt the runaway collapse of such massive clumps. These feedback effects and 
the detailed evolution of the clumps are very hard to model numerically in detail, so instead 
we adopt a more phenomenological approach, using a parametrized equation of state that 
allows us to represent a clumpy medium with the observed range of densities. The detailed 
evolution of the clumps is not relevant for our problem because the gravitational drag effect 
depends strongly on only one property of the gas, its local density. Our motivation here is 
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not to develop a fully realistic description of the gas, but only to properly model the orbital 
decay of the MBH binary caused by gravitational drag effects in a non-uniform medium. 

In order to prevent the formation of runaway clumps, we adopt a polytropic equation 
of state 

P = Kp 7 (3) 

where 7 is equal to 5/3 and K is a free parameter that we vary to produce different degrees 
of dumpiness. We consider 4 different values for K: 0.4665, 0.933, 1.3995 and 2.3325. As 
expected, the smaller K values produce a higher level of dumpiness, with denser condensa- 
tions. For K values of 0.4665, 0.933, 1.3995 and 2.3325, the gaseous disk has respectively 
60%, 30%, 20% and 10% of its total mass in regions with densities higher than 10 5 cm~ 3 . 

In the model that we just described, we introduce an MBH binary in a circular orbit, 
and we follow the subsequent dynamical evolution of the system. We evolve the system 
using the SPH code called GADGET (Springel, Yoshida & White 2001). We are interested 
in studying the evolution of the MBH binary separation with different values of the model 
parameters. The first parameter considered is the level of dumpiness of the gas, which we 
vary by varying the value of K. Secondly, since the total amount of gas in the nucleus of 
a merging system varies from one case to another, and since the black hole masses may 
also vary, we consider several different values for the mass ratio between the MBHs and the 
gaseous disk. Thirdly, since the black holes need not orbit in the plane of the disk, we vary 
the inclination angle between the plane of the binary and the plane of the disk. Table 1 lists 
the parameters used in the different simulations. 

Table 1: Run Parameters 



RUN 


M BH /M gas 


angle 


K 


A 


0.01 


0° 


0.4665 


B 


0.01 


0° 


0.933 


C 


0.01 


0° 


1.3995 


D 


0.01 


0° 


2.3325 


E 


0.01 


22.5° 


0.933 


F 


0.01 


45° 


0.933 


G 


0.01 


67.5° 


0.933 


H 


0.03 


0° 


0.933 


I 


0.05 


0° 


0.933 


J 


0.1 


0° 


0.933 


K 


0.3 


0° 


0.933 


L 


0.5 


0° 


0.933 
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4. Results 

4.1. Effects of Varying the Clumpiness 

We start by varying the level of clumpiness. We implement this by varying the constant 
K, considering 4 different values: 0.4665, 0.933, 1.3995 and 2.3325. To illustrate the effects, 
we consider first an MBH binary in which each black hole has a mass equal to 1% of the 
total gas mass, initially in a circular orbit with a binary separation of 400 pc, or half of the 
disk diameter. 

Figure 1 illustrates a representative stage in the evolution of the system for the four 
different values of K, showing the face-on density distribution in the gas disk at the same 
time in each case (t = 4.562 x 10 6 yr). The density distribution varies from relatively smooth 
spiral features to dense clumps and filaments as we reduce the parameter K. For K values 
of 0.4665, 0.933, 1.3995 and 2.3325, the disk thickness is respectively 20, 25, 32 and 40 pc. 
In our model, the disk thickness is determined by the sound speed cs = (iKp 2 / 3 ) 1 / 2 implied 
by the adiabatic equation of state used, and that ranges from 62 to 78 kms -1 for different 
values of K. In reality the sound speed cs corresponds to internal turbulent motions in the 
gas, which are observed to be of this order, as noted in section 2. Our aim is to represent the 
observed highly turbulent molecular disk (Downes & Solomon 1998) in which a turbulent 
velocity dispersion of ~ 100 km s^ 1 determines the disk thickness. 

Figure 2 shows the evolution of the binary separation in these four cases over several 
orbits. In the early evolution of the system, the separation diminishes due to the gravitational 
drag exerted by the background medium; this regime lasts until the binary separation is 
approximately 40 pc. Figure 3 shows the density distribution in the plane of the disk at time 
3.75 x 10 6 yr for run C. As was found in Paper I, the response of the gaseous medium to the 
presence of the MBHs is a lagging density enhancement and a spiral shock that propagates 
outward from each MBH. 

A notable feature of figure 2 is the marked increase of the coalescence timescale when 
the binary separation is less than 40 pc. This happens because the MBHs then come close 
enough that the wake formed by each MBH is disrupted by the gravity of the other MBH. As 
is shown in figure 4 for run C, the presence of the second MBH tends to smear out the wake, 
and therefore the gravitational drag decreases. This is not a surprising result because it is well 
known that, for both stellar and gaseous backgrounds, the dynamical friction approximation 
is no longer valid when the MBHs are at distances where the total mass of the background 
material within the orbit of the MBHs is comparable to the mass of the binary, because the 
background is then strongly perturbed by the presence of the MBHs. This physical regime 
corresponds to the transition between the gravitational drag and ellipsoidal torque regimes 



discussed in Paper I. 

The binary evolution has a strong dependence on the parameter K, and the cases with 
lower K (higher dumpiness) suffer a stronger deceleration. The explanation involves the gas 
density close to each MBH, which is higher for lower K. As a result, because the gravita- 
tional drag is stronger in a denser environment, the binary MBH separation shrinks faster 
in the cases with lower K. The effect of K on the local density is illustrated in Figure 5, 
which shows the average gas density within 20 pc around each MBH, as a function of time 
for the 4 different values of K. The figure shows that the density in the vicinity of each 
MBH varies strongly between the different cases, and this produces important differences in 
coalescence timescale, as seen in Figure 2. In run A, the drastic increase in density seen in 
Figure 5 compensates for the decrease in the efficiency of gravitational drag in the transition 
period, thereby maintaining the coalescence timescale roughly constant. But as we increase 
K (respectively in runs B, C and D), the increase of density is less, making the increase 
in coalescence timescale more obvious. In the isothermal sphere model of Paper I, the gas 
becomes increasingly concentrated around the central MBH binary, qualitatively as in run 
A, and as a result there is no distinct transition period. 

In section 5, we will show that in the later stages, the subsequent evolution of the binary 
leads to rapid coalescence in only another ~ 10 6 yr. Therefore the binary spends most of 
its time in the stages of evolution discussed above, making it important to understand the 
changes in coalescence timescale under different values of the model parameters. 



4.2. Effects of Varying the Orbital Inclination Angle 

Our second parameter is the inclination angle between the plane of the disk and the 
plane of the binary. We perform simulations with four different inclination angles: 0, 22.5, 
45 and 67.5 degrees. As in section 4.1, we consider an MBH binary initially in a circular 
orbit with MBH mass equal to the 1% of the total gaseous mass, and for these runs we adopt 
an equation of state with K=0.933, as in run B. 

Figure 6 shows the evolution of the binary separation in these four cases over several 
orbits. Although the case with i=22.5° has a coalescence timescale similar to the coplanar 
case, for inclination angles of 45 and 67.5 degrees the coalescence timescales are increased 
by factors of 3 and 4 respectively. The reason for this difference is that in these cases, the 
early evolution of the binary is driven by the gravitational drag exerted by the stellar bulge 
instead of by the gravitational drag of the gas disk. Figure 7 shows a comparison of the 
evolution of the binary separation for the case with i=45° (black curve) with the evolution 
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of the binary in the same bulge but without the gaseous disk (green curve). The agreement 
between these curves shows that the binary coalescence is driven up to time t=1.25 x 10 7 yr 
by the gravitational drag exerted by the stellar bulge. In all cases, the drag due to the gas 
disk eventually becomes dominant after the separation becomes small enough that the MBH 
binary spends more than 50% of its time in the disk, and the evolution thereafter proceeds 
much as in the coplanar case. 



4.3. Effects of Varying the Black Hole to Gas Mass Ratio 

Finally, we study the effect of varying the ratio between the mass of the MBHs and the 
mass of the gas disk. We consider four cases with black hole masses of 1%, 3%, 5% & 10% 
of the total gas mass. We introduce a MBH binary initially in a circular orbit in the plane 
of the disk, using an adiabatic equation of state with K=0.933. 

Figure 8 shows the evolution of the binary separation in these four cases over several 
orbits. In the early evolution of the system, when the separation diminishes due to the 
gravitational drag exerted by the background medium, the binary coalescence timescale 
depends on the black hole mass, with the more massive cases suffering a stronger deceleration. 
The time required for the binary separation to decrease by 75% varies from t ~4x 10 6 yr for 
Mbh = 0.01 M gas to t ~ 2 x 10 6 yr for M B h = 0.1 M gas . In units of the initial orbital period, 
this time varies from 2 to 1 orbital periods. The trend is similar to that obtained in Paper 
I for a smooth gaseous medium, and is qualitativelly consistent with the Chandrasekhar 
formula. 

The dependence of the binary coalescence timescale on the black hole mass changes when 
the binary separation arrives in the transition regime discussed in §4.1. Figure 8 shows that 
at later times there is no longer a clear dependence on the black hole mass; this is because 
as we increase the black hole mass, the wake formed is stronger, but the disruption of the 
wake by the gravity of the other MBH is also stronger, weakening the strong correlation 
with black hole mass. In fact, for the case of Mbh = 0.1 M gas , we see a small increase in the 
coalesce timescale in the transition regime. In this case the disk is starting to be globally 
perturbed by the presence of the MBHs as seen in Figure 9, and not only locally as before. 
Tidal and/or resonant forces are starting to evacuate gas from the central region, and these 
effects become more important as we continue increasing the black hole mass. As a result, 
in the overall evolution of the binary separation, the dependence on the black hole mass is 
only weak, mainly because there is not a clear dependence on mass during the later times. 

In order to determine when the gas becomes no longer important in driving the coa- 
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lescence of the binary, we continue increasing the black hole mass. We consider two more 
cases with black hole masses of 30% and 50% of the total gas mass. Figure 10 shows the 
evolution of the binary separation in these two cases over several orbits. We found that the 
coalescence tends to stall in the case where the black hole mass equals 50% of the total mass 
(red line in Figure 10). This happens because strong tidal and/or resonant forces create a 
circumbinary gap, as seen in Figure 11, and therefore the gravitational drag is no longer 
effective. The interaction of a binary with a circumbinary disk in cases where the total gas 
mass is typically equal to or smaller than the binary mass has been widely studied in the 
context of star/planetary formation (Lin & Papaloizou 1979; Goldreich & Tremaine 1980; 
Artymowicz & Lubow 1994,1996; Armitage & Natarajan 2002), where a similar circumbinary 
gap is found. The criterion for opening a gap will be discussed in §6. 



5. Final Evolution Using Higher Resolution Simulations 

The runs described above were continued until the MBH separation approaches the 
assumed gravitational softening length of 4 pc, at which point the evolution of the system 
artificially stalls. To continue the evolution of the binary it is necessary to reduce the 
gravitational softening length, something that is extremely expensive computationally. We 
choose to reduce it from e so f t =4 pc to 0.1 pc in order to follow the evolution of the binary 
separation by one more order of magnitude, and we apply this procedure to the cases where 
each MBH has 1% of the total mass in gas and the binary is in the plane of the disk, again 
for the 4 runs with different values of K. We reduce the softening length when the binary 
sepatation is approximately 30 pc, which occurs at different times for different values of K. 

Figure 12 shows the evolution of the binary separation in these four cases over several 
orbits. The black curves are the same calculations shown in Fig. 2 but on a logarithmic 
scale, and the red curves show the results for the simulations with smaller softening length. 

In the early evolution of the system, the black and red curves show almost the same 
behavior in each case. The situation changes drastically when the binary separation becomes 
less than about 10 pc, that is, when the binary arrives at separations comparable to the 
'gravitational influence' radius of the black hole: Ri n f = 2GMBH/(v BH + Cg). For example, at 
time 4 x 10 6 yr in run A, in which the black hole mass Mbh is 5 x 1O 7 M , the MBH velocity 
vbh is 172 km s -1 and the sound speed cs is (y-) 1 / 2 = 158 km s -1 in the vicinity (within 5 
pc) of each MBH, the gravitational influence radius Rj n f is 7.9 pc. At these distances the 
binary completely dominates the gravitational potential in its vicinity, and the response of 
the medium to that gravitational field is a trailing ellipsoidal density enhancement, as was 
found in Paper I. The formation of the ellipsoidal density enhancement typically occurs when 
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the binary separation is about 1.5Rj n f, and that value ranges from 9 to 12 pc in these four 
cases. Therefore the transition to the new regime is only weakly dependent on K. 

Figure 13 shows the density enhancement produced by the MBH's in the gaseous 
medium, for the simulation with K=1.3995 at four different times: 10.8, 11, 11.2 and 
11.4 x 10 6 yr. Figure 13 clearly shows that the rapid decay in the binary separation, which 
starts at around t=11.2 x 10 6 yr, coincides with the formation of the ellipsoid. In this figure, 
the black circles centered on each MBH indicate the gravitational influence radius Ri n f, and 
the figure shows that the trailing ellipsoid forms when these circles overlap. The axis of 
the ellipsoid is not coincident with the binary axis but lags behind it, and this offset pro- 
duces a torque on the binary that is now responsible for the angular momentum loss. This 
mechanism was extensively studied in Paper I and is able to reduce the binary separation to 
distances where gravitational radiation is efficient. 

In the simulations described above, we reduced the gravitational softening length but 
we maintained the number of SPH particles used. In order to verify that we have enough 
resolution to resolve adequately the region inside r = 40 pc, where the dynamics of the final 
coalescence occurs, we repeat one of these simulations increasing the number of particles. We 
use the particle splitting procedure described in Paper I to achieve this goal. The procedure 
is applied to the case where K= 1.3995. We split the SPH particles at a time of 10 7 yr when 
the binary enters the r < 30pc region; for each parent particle we introduce iV sp ii t =8 child 
particles and therefore we increase the number of particles to 1,882,648. 

Figure 14 shows the evolution of the binary separation for the high-resolution calcula- 
tion. The red curve is the same calculation shown in Fig. 12, and the green curve shows the 
MBH's separation in the high-resolution calculation. The result is qualitatively the same 
and quantitatively very similar to the low-resolution calculation. This supports the validity 
of the results shown by the red curves in Fig. 12, based on the low-resolution calculations. 



5.1. Gravitational Radiation 

The final phase in the evolution of a binary MBH occurs when eventually the MBHs 
become close enough to allow gravitational radiation to become an efficient mechanism for 
angular momentum loss. 

For an equal mass binary and scaled to convenient units, the gravitational radiation 
timescale (Peters 1964; see Eq. 16 in Paper I) is: 
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where a is the binary separation, Mbh is the mass of each black hole, and F(e) (see Pa- 
per I) contains the eccentricity dependence, which is weak for small e; e.g. F(0) = 1 and 
L(0.5)~0.205. For a binary in which the mass of each MBH is 5 x 10 7 M Q (corresponding in 
our simulations to the case where each MBH has 1% of the total mass in gas) and which has 
a separation of 0.01 pc, the gravitational radiation coalescence timescale is t gr = 1.16 x 10 7 yr. 
If we extrapolate our results to separations of O.Olpc, we predict that the MBH binary will 
merge in ~ 10 7 yr for all values of K. We feel confident that we can extrapolate our results 
to smaller separations because the ellipsoidal torque model studied in Paper I describes suc- 
cessfully the late evolution of the binary MBH, and it predicts that the separation should 
reach 0.01 pc within less than another 10 6 yr. 



6. Gap Opening Criteria 

As we found in §4.3, when the MBH binary is able to open a circumbinary gap in the 
disk, the gas no longer plays an important role in the coalescence of the binary; the problem 
then reduces to the evolution of a black hole binary at the center of a stellar system, which 
has been extensively studied in the literature (Begelman, Blandford, & Rees 1980; Makino & 
Ebisuzaki 1996; Quinlan 1996; Milosavljevic & Merritt 2001, 2003). Therefore it is important 
to determine the critical MBH mass needed to open the gap. 

The criterion for opening a gap can be determined by comparing the gap-closing and 
gap-opening times (Goldreich and Tremaine 1980). The angular momentum AL that must 
be added to the gas to open a gap of radius Ar in a disk with thickness h is AL rs p(Ar) 2 hrv. 
Using the dynamical friction formula (Chandrasekhar 1943; Ostriker 1999), we infer that the 
torque that each black hole exerts on the disk is 

(CM \ 2 
— — x /(*«.**) (M) (5) 

VBB. J 

where -M=t>BH/Cs and j( star >g as ) is a dimensionless factor which depends on the nature, 
stellar or gaseous, of the background medium (see Eqs. 7 to 10 in Paper I). This torque will 
supply angular momentum AL in a time At open = AL/T. The tendency of for gap formation 
is opposed by turbulent diffusion, which fills up a gap of width Ar on a timescale At c i osc = 
Ar/vtm-b; as mentioned before, in our adiabatic model the turbulent velocity dispersion is 
represented by the sound speed (C$ = twO- Gap formation occurs if At close > At opcn , and 
taking into account that i>bh — (GM{r)/r) 1 / 2 where M(r) is the total enclosed mass at radius 
r, this criterion can be written 



Mbe>J^i£ttti (£)^(0, ( 6 ) 



a 
4xf(M)M 
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where a is the ratio of gap radius to disk thickness, defined by Ar = ah. The simulation with 
black hole masses equal to 50% of the total gas mass suggests that a stable circumbinary gap 
forms when the MBHs are at a radius r ~ h = 25pc, with a gap width of Ar ~ 3h = 80pc 
(a — 3). For the parameters used in our simulation, Equation 6 corresponds to the condition 
Mbh > 0.46 = 0.46 M gas , in agreement with our result of Section 4.3 (Figure 11) that a 
cicumbinary gap opens when Mbh = 0.5M gas . 

As seen in Equation 6, the critical mass at a given radius depends on the total en- 
closed mass. Because approximately 80% of the dynamical mass in the region of inter- 
est is contributed by the stellar bulge (Downes and Solomon 1998), we will assume that 
M(r) = G -1 rof , a c being the central stellar velocity dispersion. The mass of the central 
black holes in nearby galaxies is observed to be correlated with the central velocity disper- 
sion, measured within the central kpc, by the so-called 'm — c c ' relation (Tremaine et al 
2002; Merritt and Ferrarese 2002): 



Or 



1 



M — = 2M - = laooi^FTj 10M - (7) 

where we assume that the mass of the binary correlates in the same way with the central 
velocity dispersion of the merged bulges. Using this relation, Equation 6 can be written in 
a simpler form: 

M bh> . 77T7VT1 (—) l-74-lO 6 M . (8) 



Anf(M)M \pc, 

As mentioned earlier, our result of Section 4.3 (Figure 11) suggests that a stable circumbinary 
gap forms when Ar ~ 3h, and therefore we set a — 3. For a Mach number of the order of 
unity, the critical mass required for a gap to form then reduces to 

/ 7 \ 2 

M BH > I— J 7.2-lO 5 M . (9) 

For typical ULIRGs (Downes and Solomon 1998), h is at least of the order of 40pc. Therefore 
the formation of a circumbinary gap will be prevented, since the observed MBH masses range 
from 10 6 to a few times 10 9 M and thus are smaller than the critical mass required to open 
a gap as given by Equation 9. The formation of a circumbinary gap will be prevented and 
the MBH binary will therefore coalesce as long as the disk thickness remains of the order of 
tens of parsecs. Observations show that even in nuclear molecular disks of mass ~ 10 8 M Q 
the velocity dispersions are still ~ 40kms _1 , indicating that disk thicknesses of the order 
of tens of parsecs are present, and this result is independent of whether or not the galaxy 
hosts a starburst (Jogee, Scotville and Kenney 2004). So gas should play a dominant role in 
driving the mergers of black holes as long as the amount of gas present in the central kpc 
is equal to or larger than this. In most galaxy mergers it is expected to have a nuclear disk 
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of at least this mass, because a disk mass of the order of 1O 8 M corresponds for an average 
spiral galaxy like the Milky Way to only 0.1% of the total mass being in gas, a small gas 
mass fraction compared to the 1 to 50% present in Sa to Scd galaxies (Young et. al. 1995). 
Taking into account that in a galaxy merger typically 60% of the gas originally present in 
the merging galaxies ends up in a massive nuclear disk, and allowing for some gas depletion 
due to star formation, we can conservately conclude that in a merger of galaxies with at 
least 1% of their total mass in gas, a nuclear disk with mass > 1O 8 M will be present, and 
that the gas will then play a major role in the evolution and final coalescence of the binary. 

On the other hand, if the masses of the black holes are greater by an order of magnitude 
or more than is expected from the 'm — <t c ' relation, it is more probable that the MBH binary 
will form a circumbinary gap inside the central lOOpc of the disk. In this case, the feeding 
of gas into the MBHs and therefore their growth will be almost stopped. This suppresses 
an important accretion phase for MBHs, in which a non-negible part of the MBH mass is 
expected to be assembled because without the gap, the accretion is predicted to be at a 
super-Eddington rate (Dopita 1997). However the bulge will continue accreting gas and 
growing, especially in the central kiloparsec where a c is measured, and therefore the stellar 
velocity dispersion will also continue growing and trying to reach the 'm — <r c ' relation. Thus 
gap formation can act as a self-regulatory mechanism for MBH growth, and this may help 
to explain the existence of the 'm — <t c ' relation. 



7. Conclusions 

In our previous paper (Paper I) we studied the role of gas in driving the evolution of a 
binary MBH, and we followed its evolution through many orbits and close to the point where 
gravitational radiation becomes important. In Paper I, we presented results for a relatively 
idealized case in which the gas was assumed to be in a nearly spherical and relatively smooth 
distribution, and we found important differences in the evolution of a binary MBH in a 
gaseous medium, compared to a stellar background. In particular, we didn't find any sign 
of ejection of the surrounding gas, as happens with stars in the later stages in the evolution 
of a binary MBH in a stellar system. 

In the present paper, we have extended this work by studying more realistic models in 
which the gas is in a disk with significant dumpiness. We do not attempt a fully realistic 
description of the gas in merging galaxies, which often exhibits starburst activity, but adopt 
a simple polytropic equation of state with a variable coefficient that allows us to represent 
gas with a varying degree of dumpiness. In this way we are able to represent gas with the 
observed range of densities in merging systems. We also vary the inclination angle between 
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the plane of the binary and the plane of the disk, and the mass ratio between the MBHs and 
the gaseous disk. 

In the early evolution of the system, as in Paper I, the separation diminishes due to 
the gravitational drag exerted by the background medium. We find a new transition regime 
that was not apparent in Paper I, when the MBHs come close enough that the wake formed 
by each MBH is disrupted by the gravity of the other MBH. In this transition regime the 
gravitational drag is less effective and therefore the coalescence timescale increases. 

In the variety of simulations that we perform, we find that varying the level of dumpiness 
of the gas, or the inclination angle between the plane of the binary and the plane of the disk, 
changes the binary coalescence timescale by less than a factor of 3. We also find that the 
dependence on the mass ratio between the MBHs and the gaseous disk is only weak, mainly 
because there is not a clear mass dependence in the transition regime. In all of the cases 
considered, the coalescence timescale varies between 5 x 10 6 yr and 2.5 x 10 7 yr, or between 
2.5 and 12 initial orbital periods. 

The final evolution of the binary, when the binary MBH dominates the gravitational 
potential in its vicinity, is driven by the same trailing ellipsoidal density enhancement that 
was found in Paper I. The axis of the ellipsoid lags behind the binary axis, and this offset 
produces a torque on the binary that is now responsible for the continuing loss of angular 
momentum. The formation of the ellipsoidal density enhancement typically occurs when the 
binary separation becomes comparable to the 'gravitational influence' radius of the black 
hole: Ri n f = 2GMbh/(v| h + Cs). As was studied in Paper I, this mechanism is able to reduce 
the binary separation to distances where gravitational radiation is efficient and will cause a 
merger within another 10 7 yr. 

We found that gas plays an important role in the coalescence of the binary when the 
MBHs satisfy the observed 'm — cr c ' relation and the disk thickness is of the order of tens of 
parsecs. This disk thickness is inferred from observations whenever gas is present in amounts 
greater than 10 8 M Q ; this corresponds, for an average spiral galaxy like the Milky Way, to 
only 0.1% of the total mass being in gas. We predict that the MBH binary will merge within 
a few times 10 7 yr for the various values of the model parameters studied. Galaxies typically 
merge in 10 8 yrs, and therefore these results imply that in a merger of galaxies with at least 
1% of their total mass in gas, the MBHs will coalesce soon after the galaxies merge. We also 
found that only if the MBHs considerably depart from the 'm — cr c ' relation is it probable 
that the MBH binary will form a circumbinary gap that stalls the coalescence, but this gap 
formation can act as a self-regulatory mechanism that helps to explain the existence of the 
'm — <t c ' relation. 
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The final coalescence has crucial implications for possible scenarios of massive black 
hole evolution and growth. Our work supports scenarios of hierarchical build-up of massive 
black holes, that in their most recent versions (Kauffmann & Haehnelt 2000; Haehnelt 2003a; 
Volonteri, Haardt & Madau 2003; Di Matteo et al. 2003) assume that the first 'seed' black 
holes appear at high redshifts (z > 10) with intermediate masses (~ 300M Q ), and that the 
black holes grow by mergers and gas accretion following the merger hierarchy from early 
times until the present. In particular, this result supports the assumption that halo mergers 
lead to black hole mergers when gas is present, something probably always true at high 
redshift. The final coalescence of the black holes leads to gravitational radiation emission 
that would be detectable up to high redshift by LISA (Menou 2003 or Haehnelt 2003b). This 
is a promising way of assessing the role of mergers in black hole growth and evolution. 
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Fig. 1. — Density distribution in the plane of the gas disk, coded on a logarithmic scale, 
at time t = 4.562 x 10 6 yr for runs A, B, C and D, which have different assumed levels of 
dumpiness. The MBHs are indicated in each plot by the black dots, and each MBH has 
1% of the mass of the gas disk. The dumpiness is varied by variying the parameter K in 



-19- 



Mgas=5 10 9 Msun Mbh = 0.01Mgas i=0 P = Kp 5 / 3 
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Fig. 2. — This plot shows the evolution of the binary separation for runs A, B, C and D 
with four different values of K, as indicated. In all cases the mass of each MBH is 1% of the 
mass of the gas, and the orbit of the binary is in the plane of the disk. A transition phase 
with a period of slower evolution after t ~5x 10 6 yr becomes increasingly prominent as K 
is increased and the gas disk becomes less clumpy. 
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Fig. 3. — Density distribution in the plane of the gas disk for run C at time t = 3.75 x 10 6 yr, 
showing in more detail the spiral shocks produced by each orbiting MBH. Qualitatively as 
was found in Paper I, the response of the gas to the presence of the MBHs is a lagging 
density enhancement and a spiral shock that propagates outward from each MBH. 
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Fig. 4. — Gas density in the plane in the transition regime, illustrated for run C at time 
t = 1.08 x 10 7 yr. The black dots indicate as before the positions of the MBHs, each with 1% 
of the mass of gas. The figure indicates how the gravitational wake of each MBH is disrupted 
and smeared out by the gravity of the other MBH, accounting for the period of relatively 
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Fig. 5. — This plot shows the average density within 20 pc around each MBH as a function 
of time for the 4 different values of K: 0.4665, 0.933, 1.3995 and 2.3325. The density in the 
vicinity of each MBH increases strongly as K decreases, accounting for the faster evolution 
shown in Figure 2 for smaller values of K. 
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Fig. 6. — This plot shows the evolution of the binary separation for runs B, E, F and G with 
inclination angles of 0°, 22.5°, 45° and 67.5° between the plane of the disk and the plane of 
the binary. For all of these runs we adopt an equation of state with K=0.933 and the mass 
of each MBH is 1% of the mass of the gas. 



-24- 



Mgas = 5 10 9 Msun K=0.933 i=45° Mbh = O.OlMgas 



400 




5 10 

time[10 6 yr] 



Fig. 7. — This plot shows the evolution of the binary separation for the case where inclination 
angle between the plane of the disk and the plane of the binary is 45° . The black curve shows 
an enlargement of the curve shown for run F in Figure 6, and the green curve shows the 
result of a run assuming the same stellar bulge but no the gas disk. 
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Fig. 8. — This plot shows the evolution of the binary separation for runs B, H, I and J 
which assume 4 different ratios between the mass of each MBH and the mass of the gas disk: 
0.01, 0.03, 0.05, 0.1. For these runs we adopt an equation of state with K=0.933, and the 
orbit of the binary is in the plane of the disk. Overall, there is only a weak dependence of 
the resulting coalescence time on the mass of the MBHs, owing mainly to the weak or even 
inverse dependence during the later stages. 
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Fig. 9. — The density distribution, in the plane of the disk for run J at t = 3 x 10 6 yr. The 
black dots indicate the position of the MBHs and in this simulation the mass of each MBH is 
equal to 10% of the mass of the gas. The gas disk is globally perturbed and forms transitory 
gaps, as a response to the presence of the MBHs. 
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Fig. 10. — This plot shows the evolution of the binary separation for runs K and L which 
assume relatively large MBH masses. The black curve shows the evolution of the binary 
when M B h = 0.5M gas , and the red curve shows the calculation with M B h = 0.3M gas . For 
these runs we adopt an equation of state with K=0.933, and the binary's orbit is in the 
plane of the disk. The evolution of the run with Mbh = 0.5M gas almost stalls at late times 
because of the dearing of a gap in the disk, as illustrated in Figure 10. 
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Fig. 11. — The density distribution, in the plane of the disk for run L at t = 3 x 10 6 yr. 
The black dots indicate the position of the MBHs and in this simulation the mass of the 
binary is equal to the mass of the gas. The gas response to the presence of the MBHs is a 
circumbinary gap created by the strong tidal and/or resonant forces near the MBHs. 
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Fig. 12. — This plot shows the final evolution of the binary separation for runs A, B, C and 
D with four different values of K. In all cases the mass of each MBH is 1% of the mass of the 
gas and the orbit of the binary is in the plane of the disk. The black curves are the same 
calculations shown in Fig. 2 but on a logarithmic scale. The red curves show the results for 
the simulations with smaller softening length. The black crosses in run C correspond to the 
four different figures shown in Fig. 13. 
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Fig. 13. — Density distribution in the plane of the gas disk for run C, illustrated at four 
different times: 10.8, 11, 11.2 and 11.4 x 10 6 yr that are indicated by crosses in Fig. 12c. 
When the regions within the gravitational influence radius of each MBH (black circles in the 
figure) merge at around 11.2 x 10 6 yr, the gas forms a trailing ellipsoidal density enhancement 
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Fig. 14. — This plot shows the final evolution of the binary separation for the case where 
K= 1.3995 (run C). The black curve is the same calculation shown in Fig. 2 but on a 
logarithmic scale. The red curve shows the result for the simulation with smaller softening 
length, the same calculation shown in Fig. 12c. The green curve shows the MBH's separation 
in the calculation that has both a smaller softening length and a much higher numerical 
resolution (N S ph =1,882,648). The result is qualitatively the same and quantitatively very 
similar to the low-resolution calculation. This supports the validity of the results shown by 
the red curves in Fig. 12, based on the low-resolution calculations. 



